home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 19 / Mac Magazin and MacEasy Magazine CD - Issue 19.iso / Wissenschaft & Technik / FFTs for RISC 1.1 / fftTest.c < prev    next >
C/C++ Source or Header  |  1996-01-30  |  3KB  |  144 lines
  1. /*  A program to test and time complex forward and inverse fast fourier transform routines    */
  2.  
  3. #include <stdio.h>
  4. #include <stdlib.h>
  5. #include <math.h>
  6. #include "fftlib.h"
  7.  
  8. #define    MACTIMING 1            /* true means time fft with macintosh calls    */
  9. #if MACTIMING
  10. #include <timer.h>
  11. #endif
  12.  
  13. #define    NUMROWS 1            /* process Matrix of NumRows different ffts of length N    */
  14. #define N 1024                /* size of FFT must be a power of 2 */
  15. #define NTIMES 20            /* number of timings,invalid if too big (if a[0][0].Re = 0|nan)*/
  16.  
  17. typedef  struct{
  18.     float Re;
  19.     float Im;
  20.     } Complex;
  21.  
  22. void main(){
  23. float        *Utbl;
  24. Complex        (*a)[N];
  25.  
  26. #if MACTIMING
  27. UnsignedWide         TheTime1;
  28. UnsignedWide         TheTime2;
  29. UnsignedWide         TheTime3;
  30. double        TheTime;
  31. #endif
  32.  
  33. long         i, il;
  34. long         TheErr;
  35. long        M;
  36.  
  37. Utbl = (float *) calloc((N/4+1),sizeof(float));
  38. if (Utbl==0)
  39.     TheErr = 2;
  40. else
  41. TheErr = FFTInit(&M, N, Utbl);
  42.  
  43. if(!TheErr){
  44.     a = (Complex (*)[N]) calloc(NUMROWS*N,sizeof(Complex));
  45.     if (a == 0) TheErr = 2;
  46. }
  47.  
  48. if(!TheErr){
  49.  
  50.             /*  set up a simple test case */
  51.     for (il=0; il<NUMROWS; il++){
  52.         for (i=0; i<N; i++){
  53.             a[il][i].Re = sqrt(il+i+.77777);    
  54.             a[il][i].Im = (il+i+.22222)*(il+i+.22222) / N - N/2;    
  55.         }
  56.         a[il][0].Re = N+3;
  57.         a[il][1].Re = 1-N;
  58.     }
  59.  
  60. #if MACTIMING
  61.     Microseconds(&TheTime1);
  62.  
  63.     for (i=0;i<NTIMES;i++){        /* do NTIMES times for timing */
  64.         ffts((float *)a, M, NUMROWS, Utbl);
  65.     }    
  66.  
  67.     Microseconds(&TheTime2);
  68.  
  69.     for (i=0;i<NTIMES;i++){        /* do NTIMES times for timing */
  70.         iffts((float *)a, M, NUMROWS, Utbl);
  71.     }    
  72.  
  73.     Microseconds(&TheTime3);
  74.  
  75.  
  76.     TheTime = (double)(TheTime2.hi - TheTime1.hi) * 65536.0 * 65536.0;
  77.     TheTime = (TheTime + (double)(TheTime2.lo - TheTime1.lo))/1000.0;
  78.     printf("Time fft = %12f  ms.", TheTime/NTIMES/NUMROWS, a[0][0].Re);
  79.     TheTime = (double)(TheTime3.hi - TheTime2.hi) * 65536.0 * 65536.0;
  80.     TheTime = (TheTime + (double)(TheTime3.lo - TheTime2.lo))/1000.0;
  81.     printf("  ifft = %12f  ms. a[0][0].Re= %6e\n", TheTime/NTIMES/NUMROWS, a[0][0].Re);
  82.  
  83. #else
  84.     printf("start timing \n");
  85.     for (i=0;i<NTIMES;i++){        /* do NTIMES times for timing */
  86.         ffts((float *)a, M, NUMROWS, Utbl);
  87.         iffts((float *)a, M, NUMROWS, Utbl);
  88.     }    
  89.     printf("end timing \n");
  90. #endif
  91.  
  92.     printf("\n");
  93.             /*  set up a simple test case */
  94.     for (il=0; il<NUMROWS; il++){
  95.         for (i=0; i<N; i++){
  96.             a[il][i].Re = sqrt(il+i+.77777);    
  97.             a[il][i].Im = (il+i+.22222)*(il+i+.22222) / N - N/2;    
  98.         }
  99.         a[il][0].Re = N+3;
  100.         a[il][1].Re = 1-N;
  101.     }
  102.  
  103.     ffts((float *)a, M, NUMROWS, Utbl);
  104.  
  105.     if (N*NUMROWS <= 256){
  106.          for (il=0; il<NUMROWS; il++){
  107.             printf("atrans = [ \n");
  108.             for (i=0; i<N; i++)
  109.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[il][i].Re, a[il][i].Im);
  110.             printf("]; \n");    
  111.             }
  112.         }
  113.     else { /* abbreviate big output */
  114.         printf("the first fft's last 32 values are: \n");
  115.         for (i=N-32; i<N; i++)
  116.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[0][i].Re, a[0][i].Im);
  117.     }
  118.  
  119.     iffts((float *)a, M, NUMROWS, Utbl);
  120.  
  121.     if (N*NUMROWS <= 256){
  122.          for (il=0; il<NUMROWS; il++){
  123.             printf("\n aitrans = [ \n");
  124.             for (i=0; i<N; i++)
  125.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[il][i].Re, a[il][i].Im);
  126.             printf("]; \n");    
  127.             }
  128.         }
  129.     else { /* abbreviate big output */
  130.         printf("\n the first ifft's last 32 values are: \n");
  131.         for (i=N-32; i<N; i++)
  132.                     printf(" %+20.15e + j * ( %+20.15e ) \n", a[0][i].Re, a[0][i].Im);
  133.     }
  134.  
  135.     free (a);
  136.     free (Utbl);
  137.     return;
  138. }
  139. else
  140. if(TheErr==2)    printf(" out of memory ");
  141. else    printf(" error ");
  142. return;
  143. }
  144.